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ABSTRACT 

V 

t 

v  }> 

'The  method  of  characteristics  is  applied  to  the  set  of  equations  which 
governs  the  propagation  of  axially  symmetric  torsional  shear  waves  in  non- 
homogeneous  elastic  media.  The  wave  velocities,  characteristic  equations, 
and  the  equations  governing  the  propagation  of  abrupt  changes  (discontinuous 
wave  fronts)  are  derived  in  closed  form.  Numerical  integration  along  the 
characteristic  directions  was  carried  out  for  several  examples  on  an 
electronic  computer.  The  solutions  of  three  specific  examples  calculated 
show  general  agreement  with  existing  solutions  by  other  methods.  For  certain 
problems,  the  method  of  characteristics  yield  additional  results  which  cannot 
be  obtained  by  the  Laplace  transform  method. ^  , 


NOMENCLATURE 

c  =  (G/p)1/2  *  shear  (or  distorsional)  wave  velocity 

G  =  shear  modulus  (function  of  i 
G  =  G/G 

o 

r  =  radial  distance 
rQ  =  inner  radius  of  plate 
r  =  r/r 

o 

t  *  time 

t  =  c  t/r 

o  o 
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v  =  tangential  displacement 

vt  =  3v/3t  *  particle  velocity 

v  »  3v/3r  «  shear  strain 

r 

v  =  Gv/rr 

O  0  0 

p  a  density  of  material  (function  or  r) 

P  *  p/pq 

t  »  torsional  shear  stress 

?  -  T/To 

Subscript 

o  =  properties  at  rQ 


INTRODUCTION 

The  problem  of  shear  wave  propagation  due  to  a  suddenly  applied  rotary 
disturbance  in  a  homogeneous  elastic  plate  was  solved  by  Goodier  and  Jahsman1, 
The  corresponding  problem  in  a  nonhomogeneous  plate  was  solved  by  Sternberg 
and  Chakravorty2 .  In  both  these  cases  the  Laplace  transform  technique  was 
used.  Except  in  a  few  cases  of  special  radial  distributions  of  the  shear 
modulus,  all  their  solutions  are  in  integral  form  which  can  be  evaluated  only 
by  numerical  integration.  A  solution  obtained  by  the  Laplace  transform 
technique  is  applicable  for  only  one  type  of  initial  and  boundary  conditions. 
To  solve  for  a  different  type  of  condition,  the  problem  must  be  reinitiated 
and  techniques  for  inversion  developed.  In  their  study  of  the  nonhomogeneous 
plate  problems  Sternberg  and  Chakravorty  were  mainly  interested  in  the  quali¬ 
tative  effect  of  the  variation  of  shear  modulus;  therefore,  a  simple 
exponential  variation  was  selected.  Solutions  for  other  types  of  radial 
distribution  of  the  shear  modulus  are  not  available.  For  these  reasons  there 
is  a  need  for  other  methods  in  treating  shear  waves  in  nonhomogeneous  plates. 
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In  this  paper  the  propagation  of  cylindrical  shear  waves  in  nonhomogeneous 
elastic  bodies  is  treated  by  the  method  of  characteristics.  By  using  this 
method,  the  distribution  of  wave  velocity  (physical  characteristic),  the 
characteristic  equations,  as  well  as  the  equations  governing  the  propagation 
of  abrupt  change  in  stress  (step  input) ,  may  be  determined  in  closed  form. 
Numerical  integration  for  the  determination  of  the  stress  field  behind  the 
wave  front  may  be  accomplished  readily  for  any  type  of  input;  and  for  any 
type  of  radial  distribution  of  the  shear  modulus  and  density.  As  examples, 
materials  with  simple  exponential  distribution  of  the  shear  modulus  under 
step  input  in  stress  are  presented.  For  a  certain  class  media,  the  Laplace 
transform  method  yields  results  only  up  to  a  certain  critical  time;  whereas 
the  method  of  characteristics  yield  solutions  beyond  this  critical  time. 

In  Ref.  3,  the  method  of  characteristics  was  applied  to  cylindrical  and 
spherical  dilatational  waves  in  a  homogeneous  elastic  material.  The  present 
paper  is  an  extension  of  the  method,  not  only  to  the  case  of  shear  waves,  but 
also  to  nonhomogeneous  media. 


GOVERNING  EQUATIONS 

The  governing  equations  in  cylindrical  coordinates  for  elastic  torsional 
shear  waves  under  axisymmetrical  loading  conditions  are, 


3r  ^  2t  32v 
Sr  *  ~  "  p  9tT 


,3v  V, 

T  *  G  ( 3?  ~  7> 


(1) 


(2) 


where  r  is  the  radius;  t  time;  p  is  the  density;  G  is  the  shear  modulus; 

t  is  the  torsional  shear  stress;  and  v  is  the  tangential  displacement.  The 
shear  modulus  and  the  density  are  in  general  arbitrary  functions  of  the  radius. 
Substituting  eq,  (2)  into  eq,  (1),  we  obtain  a  single  second  order  equation 


•3- 


S2V  +  irVv  +  1 r 3V  _  V , 
It?  3rV  G  3r  l3r  rJ 


1  32v 

critr 


(3) 


where  c  *  (G/p) 1 /2  is  the  torsional  shear  wave  velocity.  It  may  be  noted 
that  since  only  shear  stress  is  involved,  these  equations  are  exact  for  both 
a  plate  (plane  stress)  and  a  hollowed  infinite  body  (plane  strain) . 

In  the  application  of  the  method  of  characteristics,  we  may  use  either 
eqs.  (1)  and  (2) t  the  stress  approach;  or  eq,  (3),  the  displacement  approach. 
The  governing  equations  for  both  approaches  will  be  given  below;  while  the 
numerical  procedures  for  the  displacement  approach  only  will  be  presented. 


CHARACTERISTIC  EQUATIONS 

In  the  displacement  approach,  we  apply  the  method  of  characteristics  to 
the  single  second  order  equation,  (3),  and  obtain  the  following  two  physical 
characteristics 

£  -  *  «W1/2  (4) 

which  will  be  called  the  I+  and  i"  characteristics,  respectively.  Notice  that 
eqs,  (4)  are  of  the  same  form  for  both  homogeneous  and  nonhomogeneous  materials. 
For  homogeneous  materials,  the  physical  characteristics  are  two  families  of 
straight  lines  of  constant  slope,  whereas  for  nonhomogeneous  materials,  they 
are  two  families  of  curved  lines  in  the  r,t-plane.  In  both  cases,  once  the 
distribution  of  G  and  p  are  given,  the  physical  characteristics  are  determined 
independent  of  the  loading  and  solution  of  the  problem. 

The  characteristic  equations  of  (3),  with  v  for  3v/3t,  v  for  3v/3r,  are, 

t  X 

d(vt)  +  c  d(vr)  »  t  c(vr  -  £)  (^  +  ^|)  (5) 

along  the  I*  and  i”  characteristics,  respectively.  For  homogeneous  material, 
the  term  containing  dG  vanishes  in  (5) , 

In  the  stress  approach,  we  differentiate  eq.  (2)  with  respect  to  time  and 
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rewrite  (1)  and  (2)  as 


3t 

3r 


+ 


3x 

3t 


3v  v 
_  /  t 

G  c~  ■  r} 


(6) 

(7) 


These  may  be  considered  as  two  first  order  equations  in  terms  of  r  and  v^.. 
Applying  the  method  of  characteristics  to  eqs.  (6)  and  (7),  we  obtain  two 
physical  characteristics  identical  to  (4) .  The  corresponding  characteristics 
equations  are 

dt  +  —  d(v  )  =  (-2t  +  —  v  )  —  (8) 

c  t'  c  t  r  ' 

along  I+  and  I" ,  respectively.  Unlike  eqs,  (5),  the  form  of  eqs,  (8)  is  the 
same  for  both  homogeneous  and  nonhomogeneous  materials.  It  can  be  shown 
readily  that  upon  substitution  of  (2)  into  (8),  eqs.  (5)  are  obtained. 

For  the  present  problem,  the  stress  approach  and  the  displacement  approach 
yield  identical  results.  This  is  different  from  the  problems  of  spherical  and 
cylindrical  dilatational  waves3  and  the  problem  of  cylindrical  flexural  waves 
in  a  plate4.  In  both  those  cases,  the  stress  approach  produces  one  extra 
physical  characteristic,  dr/dt  »  0,  which  has  an  associated  characteristic 
equation  equivalent  to  a  restatement  of  the  static  stress-displacement  relations. 


PROPAGATION  OF  DISCONTINUITY 

Across  the  physical  characteristics  the  second  derivatives  of  v  (or  the 
first  derivatives  of  t  and  v  )  may  be  discontinuous.  Discontinuities  of  the 
first  derivatives  of  v  (or  t  and  vt  themselves)  may  also  exist  across  the 
physical  characteristics,  but  these  will  not  be  governed  by  eqs.  (5)  or  (8). 

In  Ref.  3,  the  equations  governing  the  discontinuities  in  the  first  derivatives 
(jump  conditions)  of  the  displacement  variable  in  dilatational  waves  are 
derived  by  using  the  stress  approach.  In  Ref,  4,  a  similar  set  of  jump 
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conditions  are  derived  for  flexural  waves  by  using  the  displacement  approach. 
In  this  paper,  we  shall  follow  the  displacement  approach  and  derive  the  jump 
conditions  for  shear  waves  in  nonhomogeneous  materials. 


Let  A  and  B  be  two  points  on  a  I  characteristic  as  shown  in  Fig,  1. 

The  two  I*  characteristics  passing  through  A  and  B  are  represented  by  I *  and 
I21  respectively.  If  a  discontinuity  of  vt  across  l|  exists,  then  v  A  -  vtJJ  ** 
is  finite  but  different  from  zero  as  l!|  is  allowed  to  approach  I*,  or  as  dr 
approaches  zero.  Writing  eq.  (5)  (with  the  lower  sign,  for  I”)  and  integrating 
from  A  to  B,  we  have 


+ 


fB 

c  d(v  ) 
'A  r 


v,dr 
rJ  r 


rB 

c(v 

JA 


v.dG 
rJ  G 


(9) 


As  B  approaches  A,  the  right  hand  side  of  (9)  vanishes,  since  the  integrands 
contain  bounded  values  of  c,  v  ,  and  v,  provided  the  domain  in  the  physical 
plane  does  not  include  the  line  r  =  0.  Integrating  the  left  hand  side  integral 
by  parts,  and  keeping  in  mind  that  vr  is  bounded,  as  B  approaches  A  eq.  (9) 
becomes 

[vj  +  cEvr]  *  0  do) 


where  brackets  are  used  to  designate  jumps.  The  variations  of  [v^]  and  [vy] 
as  they  propagate  along  the  I  is  obtained  by  writing  (5) ,  with  the  upper 
signs,  along  I*  and  I*,  and  subtracting  one  from  the  other.  As  B  approaches 
A,  we  have 

d[vt:i  -  C  d[vr]  -  c[vrl  &  *  c[vr]  f  (11) 

where  the  condition  [v]  «  0  has  been  utilized.  Eliminating  [v  ]  from  (11)  by 
(10),  we  obtain 


-  7  c— ♦ 

2  v  r 


(12) 


which  may  be  integrated  to  give 


Equations  (10)  and  (13)  then  yield 

[vr]  "  +K  C^1/2  (14> 

The  corresponding  jump  in  t  obtained  from  (2)  and  (14)  is  governed  by 

H  *  >/2  (is) 

Equations  (13)  to  (15)  are  for  jumps  across,  and  the  propagation  along, 
a  I+  characteristic.  Following  the  same  procedure,  equations  for  those  across 
and  along  a  1“  characteristic  can  be  shown  to  be 

lyt]  -  «  (^)i/2 

Cvr]  *  (gHt>  1/2  C16> 

W  -  *K(S-)l/2 

The  same  set  of  equations  (13)  to  (16)  may  also  be  derived  from  the  stress 
approach , 


INITIAL  AND  BOUNDARY  CONDITIONS 

The  elastic  body  under  consideration  will  be  either  an  infinite  sheet 

with  a  circular  hole,  or  an  infinite  hollow  cylinder.  These  configurations 

can  be  represented  by  rQ  <  r  <  «,  where  rQ  is  a  constant.  Initially,  the  body 

is  not  loaded,  thus  the  stress  and  velocity  are  zero.  For  time  greater  than 

zero,  the  input  is  applied  at  the  boundary  r  *  rQ,  either  suddenly  or 

gradually.  This  input  can  be  in  the  form  of  specified  time  functions  of  any 

one  of  the  three  variables,  v„,  v  ,  or  t, 

t  r 

NUMERICAL  PROCEDURE 

It  is  convenient  to  introduce  non-dimensional  quantities  as  follows: 
i  =  r/r0,  t  .  tco/ro,  5  •  G/G0,  7  =  V/Toro>  *  *  t/V 

(17) 

p  ■  p/pQ,  c  *  c/cQ 
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Thus,  the  results  in  terms  of  these  quantities  are  true  for  materials  of 

any  values  of  G  ,  p  ,  and  r  . 

o*  o*  o 

A  numerical  technique  for  stepwise  integration  along  the  physical 
characteristics  is  developed.  In  the  r,t-plane,  the  region  between  r  =  1 
and  r  *  1  +  ct  is  divided  into  a  grid  system  by  the  two  families  of 
characteristics.  The  displacement  approach  is  used,  thus  at  each  grid 
point  values  of  v,  vt,  and  vr  are  calculated.  Since  only  continuous  v  is 
considered  for  all  regions  in  the  physical  plane,  we  may  write  the  continuity 
equation 

dv  =  v.  dt  +  v  dr  (18) 

t  r 

along  any  directions.  In  our  numerical  work,  this  continuity  equation  along 

the  I"  characteristics  is  used.  Values  of  the  three  variables  v,  v,.  and  v 

t  r 

at  a  typical  interior  point  1  of  Fig.  2  may  be  calculated  from  eqs,  (5)  and 
(18)  expressed  in  finite-difference  form,  if  all  quantities  at  the  neighboring 
points  2  and  3  are  known.  Along  the  leading  I+  characteristic  passing  through 
(1,0),  all  three  variables  vanish  if  the  prescribed  boundary  condition  at  (1,0) 
is  continuous.  For  jump  input  at  (1,0)  values  of  v  and  v  along  the  leading 
characteristic  are  calculated  from  eqs.  (13)  and  (14).  Along  the  boundary 
r  =  1,  either  v  or  t  is  specified;  correspondingly,  the  I+  characteristics 

I* 

to  the  left  are  absent  leaving  two  equations  for  two  unknowns. 

In  the  numerical  calculation,  the  characteristic  grid  system  was  constructed 
by  choosing  points  on  the  leading  I+  characteristic  with  equal  horizontal 
distance,  as  shown  in  Fig,  2,  The  I"  characteristics  are  constructed  from 
the  reflections  of  the  I+  characteristics  from  the  boundary  r  »  1, 


SPECIFIC  EXAMPLES 

A  few  specific  examples  of  various  inputs  at  r  =  1  are  calculated  and  the 
results  compared  with  existing  solutions  by  other  methods.  Although  all  the 
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equations  and  the  numerical  procedures  discussed  in  this  paper  are  applicable 
to  bodies  with  arbitrary  radial  distribution  of  6  and  p,  in  the  specific 
examples  presented  below  a  special  distribution  is  selected,  i.e,, 

G  *  ra ,  p  *  1 

(19) 

3  .  CG/p)1/2  «  ra/2 

where  a  is  a  constant.  With  these  functions  of  G  and  p,  our  results  may  be 
compared  directly  with  those  of  Refs.  1  and  2.  For  a  unit  step  stress  input, 
the  long  time  asymptotic  solution  of  stress  should  approach  the  corresponding 
static  solution,  i.e.,  x  *  l/rz,  regardless  of  the  elastic  properties  of  the 
medium. 

Homogeneous  Medium 

For  a  homogeneous  medium,  ot  *  0,  results  of  our  calculation  are  shown  in 
Fig.  3,  for  a  unit  step  x  input  at  the  hole.  On  this  figure,  the  results 
obtained  by  Goodier  and  Jahsman1  is  also  shown.  Figure  3  is  a  plot  of  x 
against  time,  at  three  different  radial  locations.  As  far  as  the  arrival 
time,  the  magnitude  of  the  peak  stress,  and  the  asymptotic  static  values  of 
the  stress  are  concerned,  our  results  and  those  of  Ref,  1  are  in  agreement. 
However,  a  slight  discrepancy  in  x  exists  during  a  time  period  after  the 
arrival  of  the  wave  front. 

It  is  interesting  to  note  that  for  homogeneous  media,  the  governing 
equation  in  terms  of  displacement,  eq.  (3),  is  of  the  same  form  as  the 
corresponding  equation  for  cylindrical  dilatational  waves,  eq,  (10)  of  Ref.  3. 
If  the  input  at  r  *  1  is  in  terms  of  prescribed  velocity,  then  the  solutions 
(displacement  and  velocity)  for  the  dilatational  wave  can  be  used  as  those 
for  the  shear  wave,  if  the  value  of  the  wave  velocity  is  properly  adjusted. 

The  stresses  must  be  calculated  from  the  proper  stress-displacement  equations 
for  each  case  separately. 
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Nonhomogeneous  Medium,  g  «  1 


With  a  value  of  a  ■  1,  the  physical  characteristics  are  curved  lines  as 
shown  in  Fig.  2.  With  a  constant  increment  in  r  along  the  leading  I+ 
characteristic,  the  grid  segments  along  the  i”  characteristics  are  not  of 
constant  length.  The  change  in  length  of  the  segments  is  not  severe,  and 
is  believed  to  be  tolerable  in  the  numerical  calculation  within  the  region 
of  interest.  Our  calculated  stress  distribution  is  shown  in  Fig.  4,  On 
this  figure,  the  results  obtained  by  Sternberg  and  Chakravorty2  are  also 
given  for  comparison.  Again,  correct  values  of  arrival  time,  peak  stress, 
and  long  time  asymptotic  stress  are  obtained  by  both  methods.  A  slight 
discrepancy  exists  for  stress  during  a  time  period  after  the  arrival  of  the 
wave  front. 

Nonhomogeneous  Medium,  a  «  10 

This  is  a  case  of  particular  interest  because  an  exact  closed  form 
solution  exists.2  For  a  unit  step  stress  input,  closed  form  solutions 
exist  for  those  media  with  the  following  values  of  a, 

a  =  24^  (k  -  0,  ±1,  ±2,  ...)  (20) 

where  a  »  10  is  a  special  case  corresponding  to  k  *  -2.  As  shown  in  Ref,  2, 

the  Laplace  transform  of  the  displacement  is  a  modified  Bessel  function  which 

degenerates  into  elementary  functions  if  a  is  given  by  (20).  These  elementary 

functions  can  be  inverted  into  closed  form  functions  of  r  and  t. 

For  problems  with  a  >  2,  the  Laplace  transform  method  yields  a  critical 

time  t  ,  where  solutions  exist  only  in  the  time  interval 
00  9 

0  <  l  <  ta  (21) 

We  shall  show  that  solutions  above  this  critical  time  can  be  obtained  from 
the  method  of  characteristics. 
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Integrating  the  dimensionless  physical  characteristic  with  c  given 
by  (19),  we  have  the  equation  for  the  leading  I+  characteristic 

t  =  {?(2“°°/2  -  1}  (22) 

This  characteristic  has  an  asymptote  at  t  **  t  »  2/  fa -2) »  where  r  *►  ®. 

In  Fig.  5,  this  leading  I+  characteristic  for  a  =  10  is  labeled  as  curve  OA. 

As  r  ->  «,  the  wave  speed  c  and  the  stress  r  all  approach  infinity.  In  the 
Laplace  transform  method,  the  transformed  displacement  is  governed  by  an 
ordinary  differential  equation  with  r  as  the  independent  variable,  which 
cannot  tolerate  unbounded  boundary  values.  Consequently,  no  solution  can  be 
obtained  for  t  >  t  from  the  Laplace  transform  method. 

From  the  principle  of  domain  of  dependence  in  the  method  of  character¬ 
istics,  if  proper  values  of  v  is  prescribed  along  OB  (Fig.  5),  which  is  not 
a  characteristic,  and  proper  values  of  v  and  v^  are  prescribed  on  OA,  then 
the  solution  is  uniquely  determined  in  the  region  OAB,  where  AB  is  the  I" 
characteristic  passing  through  A  and  B,  Notice  that  point  B  is  at  a  time 
larger  than  t^,  which  is  0,25  for  a  =  10.  Along  the  leading  I+  characteristic, 

values  of  v  and  v  increase  without  bound  as  r  increases.  The  domain  with 
r 

which  the  solution  can  be  obtained  is  therefore  bounded  by  line  CD,  the  I" 
characteristic  asymptotic  to  the  line  t  =  t^.  In  applying  the  numerical 
integration  along  characteristics,  accurate  solution  cannot  be  obtained  for 
points  very  close  to  the  line  CD;  because  the  characteristic  grids  are  greatly 
distorted  for  large  r  and  the  values  of  v  and  vr  are  too  large. 

Results  of  our  calculation  and  those  of  Laplace  transform  are  shown  in 
Fig.  6,  in  the  form  of  t  agains'1  t  at  different  radii.  For  t  <  t  ,  the  Laplace 
transform  solution  indicates  t  is  constant  for  fixed  radius.  The  results  from 
method  of  characteristics  are  in  complete  agreement  with  this  for  r  «*  1.1  and 
1.2.  The  Laplace  transform  solution  stops  abruptly  at  t  *  0.25;  whereas  the 
method  of  characteristics  yields  results  beyond  this  time,  and  the  stress  r 
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remains  constant  for  t  >  0,25  at  r  *  1.1,  For  large  values  of  r,  or  for 
values  of  t  close  to  the  curve  CD,  our  results  show  an  increase  in  stress, 
which  is  probably  due  to  the  inaccuracy  introduced  by  the  large  distorted  grids. 

DISCUSSION 

For  the  case  of  a  *  0  and  a  *  1,  the  numerical  results  from  the  method 
of  characteristics  deviate  slightly  from  the  curves  presented  in  Refs.  1  and 
2\  which  were  obtained  by  the  Laplace  transform  technique.  The  curves  presented 
in  Refs.  1  and  2  were  in  rather  small  scale  without  enough  resolution  for 
accurate  evaluation.  Therefore,  the  discrepancy  may  be  due  to  either  inaccuracy 
in  plotting  and  reading  of  curves;  or  inaccuracy  in  one  or  both  of  the  numerical 
results.  The  basic  procedure  of  the  present  calculation  remains  unchanged  for 
different  values  of  a;  and  for  the  case  of  a  =  10  the  present  calculation  is 
very  accurate.  This  seems  to  give  a  certain  degree  of  confidence  in  the 
accuracy  of  the  present  method. 
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